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ABSTRACT 

We report on the results of three convective dynamo simulations with an outer coronal layer. The mag- 
netic field is self-consistently generated by the convective motions beneath the surface. Above the convection 
zone we include a polytropic layer that extends to L6 solar radii. The temperature increases in this region 
to 8 times the value at the surface, corresponding to ^ 1.2 times the value at the bottom of the spherical 
shell. We associate this region with the solar corona. We find a solar-like differential rotation with radial 
contours of constant rotation rate, together with a solar-like meridional circulation and a near-surface shear 
layer. This spoke-like rotation profile is caused by a non-zero latitudinal entropy gradient which violates the 
Taylor-Proudman balance via the baroclinic term. The lower density stratification compared with the Sun leads 
to an equatorward return flow above the surface. The mean magnetic field is in most of the cases oscillatory 
with equatorward migration in one case. In other cases the equatorward migration is overlaid by stationary or 
even poleward migrating mean fields. 

Subject headings: Magnetohydrodynamics - convection - turbulence - Sun: dynamo - Sun: rotation - Sun: 
activity 



1. INTRODUCTION 

The Sun has an activity cycle of eleven years, which is man- 
ifested by sunspots occurring at the solar surface. The sunspot 
number changes from a few during minimum to over 200 dur- 
ing maximum. The sunspot locations display a latitudinal 
dependence during the cycle. At solar minimum, sunspots 
emerge preferably at higher latitudes and during maximum 
at lower latitudes. By plotting these for several cycles, one 
obtains the 'butterfly diagram' . Every eleven years the polar- 
ity of the sunspot pairs changes sign, manifesting the 22-year 
magnetic cycle. To understand this cyclic behavior, one has 
to connect the fluid motions in the Sun with magnetic field 
generation to construct dynamo models. These dynamo mod- 
els should be able to reproduce the 22-year magnetic activity 
cycle as well as the large-scale magnetic field evolution at the 
surface of the Sun. It is widely believed that the sunspots 
are correlated with the large-scale magnetic field distribution. 
Therefore, a successful solar dynamo model should reproduce 
the equatorward migration of the large-scale field as we ob- 
serve it indirectly from sunspots and more directly from syn- 
optic magnetograms. 

Until recently, only kinematic mean-field models, where 
turbulent effects are parameterized through transport coeffi- 
cients (e.g. Krause & Radler 1980), have been able to show 
this (e.g. Dik pati & Charbonneau 1999; Kapyla et al. 2006; 
iKitchatinov & Qlemskoy 20121). These models have been 
used to reproduce certain feature s of the solar cycle, such 
as the Maunder minimum (e.g., iKarakI 120101) . However, 
these models are only valid in the kinematic regime where 
the fluid motions are assumed given, so they are not self- 
consistently generated. The backreaction from the magnetic 
field is either ignored or taken into account in a rudimen- 
tary way involving ad hoc quenching for the turbulent trans- 
Electronic address: [j oem@nordita.org (Revision: 1.175 ) | 



port coefficients. Direct numerical simulations (DNS) of the 
solar dynamo have been unsuccessful in producing equa- 
torward migration using convective motio ns to dr ive a dy- 
namo (e.g.LG ilman 1983; Brun et al. 2004 lKapvla et aklMlOl: 
iGhizaru et al. 2010: Brown et al. 2011) so far. This was pos- 
sibly due to the low fluid and magnetic Reynolds numbers in 
the simulations. Equatorward migration of the mean magn etic 
field was, for the first time, found in DNS by Kapyla e t al.l 
(2012). The exact cause is not yet fully understood, but the 
amo unt of density stratification seems to play an important 
role (iKapvla et al.il2013h . 

An important ingredient of the solar dynamo is differen- 
tial rotation. It is believ ed that strong shear a t the bottom 
of the convection zone ( Spiegel & Weiss 1980) or near the 
surface (Brandenburg 2005), plays an important role am- 
plifying the magnetic field. However, even today it is not 
straightforward to reproduce a sol ar-like differential rota- 
tion profile. Mean-field simulations (iBrandenburg et al.ll 19921: 
IKitchatinov & Rtidiger T995) have been able to reproduce a 
solar-like rotation profile by modeling small-scale e ffects by 
mean-field coefficients through the A effect (see e.g. lRiidigen 
■ 1980. 1989). These models reproduce the positive (nega- 
tive) latitudinal gradient of angular velocity in the northern 
(southern) hemisphere - i.e., the equator rotates faster than the 
poles - together with 'spoke-like' contours in the meridional 
plane. DNS of convective dynamos are able to reproduce a 
rapidly rotating equator at sufficiently large Coriolis numbers 
dBrun et al.i 2004; K apvla et al] 12011 bl) . To generate spoke- 
like differential rotation in DNS, however, has so far been re- 
ported to be possible only by imposing a latitu dinal entropy 
gradi ent at the bottom of the convection zone ( Miesch et al.l 
|2006|) . A self-consistently generated spoke-like profile has 
not yet been found. 

An important issue with solar dynamo models is the ef- 
fect of catastrophic quenching of the dynamo at high mag- 
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netic Reynolds numbers; see iBrandenburg & SubramanianI 
(l2005h . This is caused by the accumulation of magnetic he- 
licity in the dynamo region. DNS provide evidence that mag- 
netic helicity fluxes both within and through the boundaries 
of the dynamo domain can prevent the dynamo from being 
catastrophically quenched (e.R . Brandenburg & Sandin 2004; 
iHubbard & Brandenburgll2()12|) . In the case of the Sun, mag- 
netic helicity flux can emerge through the solar surface and 
can be transported away fr om the Sun by coronal mas s ejec- 
tions or by the solar wind (Blackma n & Brande nburg "2003). 
In earlier work we mimic the case of the Sun by using an up- 
per layer on the top of a dynamo region to allow for ma gnetic 
helicity fluxes leaving the domain (Wamecke & BrandenburgI 
|20TqI: IWamecke et alJlIoTTl IMIallbh . This two-laver model 
was successful in showing that the dynamo is not only en- 
hanced, but it can actually trigger the emergence of coronal 
ejections. These ejections have a similar shape as coronal 
mass ejections and carry a significant amount of magnetic he- 
licity out of the dynamo region. 

In this work we use the two-layer approach to investigate 
the influence of the coronal layer as an upper boundary condi- 
tion for a convective dynamo. We focus on the physical prop- 
erties and dynamics in the convection zone. The effects of 
varying the strength of stratification on a conve ctive dynamo 
witho ut corona is studied in a companion paper (|Kapyla et al.l 
l20T3h . 

2. MODEL AND SETUP 

We use a two-layer model in spherical polar coordinates 
(r, 6>, (/)), where the lower layer (r < R) represents the con- 
vection zone and the upper one the corona. The simula- 
tion is performed in a spherical wedge with radial extent 
^0 ^ r ^ = 1.6 i?, where ro = 0.7 R corresponds 
to the bottom of the convection zone and R to the solar ra- 
dius, for colatitudes 15° < < 165° and an azimuthal extent 
< (/) < 45°. We solve the following equations of compress- 
ible magnetohydrodynamics. 
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where the magnetic field is given by ^ = V x A and thus 
obeys V B = at afl times, J = ijlq^V x B is the 
current density, /io is the vacuum permeability, r] and u are 
the magnetic diffusivity and kinematic viscosity, respectively, 
D/Dt = + n • V is the advective time derivative, p is 
the density, and u is the velocity. The traceless rate-of- strain 
tensor is given by 



(5) 



where sem i colons denote co variant differentiation; see 
iMitra et al.l (l2009h for details. Furthermore, Qq = 
1^0 (cos 6>, — sin6>, 0) is the rotation vector and p is the pres- 
sure. The gravitational acceleration is given by 



g = -GMr/r^, 



where G is Newton's gravitational constant and M is the mass 
of the star. The radiative and subgrid scale (SGS) heat fluxes 
are defined as 



prad 



7.SGS 



XsGSpTVs, 



(7) 



where K is the radiative heat conductivity and XSGS turbulent 
heat conductivity, which represents the unresolved convective 
transport of heat. The fluid obeys the ideal gas law, p = (7 — 
l)pe, where 7 = cp/cy = 5/3 is the ratio of specific heats at 
constant pressure and constant volume, respectively, and e = 
cyT is the internal energy density, defining the temperature T. 
Finally, Fcooi is the radial cooling profile which is specified 
in Equation (flQI) . 

The two-layer model is similar to that used in the pre- 
vious work (Wa mecke & Brandenburgll2010l IWamecke et al.l 
201 IL l2012atibi) . except that here we improve the model of 



Warnecke et al.l (|2012b|) in two important ways. First, we 



use a more r ealistic mode l for the convection zone as in 
iKapvla et al.l IMTaL IMl) . Instead of using a polytropic 
setup with m = 1, we lower the radiative flux by decreas- 
ing K and introducing a turbulent hea t conductivity xsGS, 
(referred to as xt in Ka pyla et al.ll2012|) . We apply a piece- 
wise constant profile for xsGS such that in the interval of 
0.75 R < r < 0.97 R it is equal to a quantity XsGS (whose 
value is similar to that of the microphysical viscosity v), 
and it goes smoothly to zero above and below the bound- 
aries of the interval. Additionally, we change the temperature 
profile compared to our earlier isothermal cold corona to a 
temperature- stratified corona, which is ^ 8 times hotter than 
the surface and ^ 1.2 times than bottom of the convection 
zone. In Figure [T] the averaged temperature profile is shown 
together with the averaged density, pressure and entropy pro- 
files for one of the runs. 

We initialize the simulations with pre-calculated radial pro- 
files of temperature, density, and pressure. In the convection 
zone (r < R) we have an isentropic and hydrostatic initial 
state for the temperature, whose gradient is given by 



dT 
dr 
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cv(7 - l)(m + 1) ' 



(8) 



where m = 1.5 is the polytropic index. This leads to a tem- 
perature minimum Tmin above the surface of the convective 
layer air = R. In the corona (R < r < Rc),we prescribe the 
temperature as 



T(r) = 



1 + tanh 



rtr 



w 



(9) 

where Tcor is temperature in the corona, and rtra and the 
width w = 0.02 are chosen to have a smooth tempera- 
ture profile as shown in Figure [T] The radial cooling profile 
TcooK^) in Equation maintains the temperature profile. 



- cool 



Tof{r){cl-c%)/cl,, 



(10) 
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where /(r) is a profile function equal to unity in r > R and 
going smoothly to zero in r < R, and Fq is a cooling lumi- 
nosity chosen so that the sound speed in the corona relaxes 
towards the temperature profile given in Equation (|9}. The 
stratification of density follows from the hydrostatic equilib- 
rium. The radial heat conductivity profile is chosen so that 
the energy in the convection zone is mostly transported by 
convective motion. We apply a profile for the viscosity u 



3 



2 







-1 





/' 

/' 




; 


P/Po 

- : 


; / 

1 

il 




; 

: 


- / ■ / 
: ^/S : / 
i /' 


/ 

/' 0.100 
0.010 
0.001 


^"^"^^r^^^ : At/To 

- ■ • , ^ ■ P/Po 
p/Po \ 








0.8 1.0 1.2 1.4 1 

, 1 , , ^/^ 1 , , , 


61 


0.8 1.0 




1.2 1.4 


1 



.6 



r/R 



Fig. 1. — Radial profiles of stratification for Run A. The normalized 
density p/po (dashed lines), pressure p/po (dash-triple-dotted), tempera- 
ture T/To (solid) are plotted together with the specific entropy s/cp (dash- 
dotted) over the radius. The inset shows various profiles in logarithmic rep- 
resentation to emphasize the steep decrease of the pressure and density in the 
coronal layer. The index represents the value at the bottom of the convec- 
tion zone. 



which is constant in the convection zone (r < R) and in- 
crease smoothly above the surface to a value that is 20 times 
higher in the corona. This helps to suppress high velocities 
and sharp flow structures aligned with the rotation vector in 
the corona especially in the beginning of the simulation, when 
the magnetic field is weak. We initialize the magnetic field 
with weak Gaussian-distributed perturbations inside the con- 
vection zone. 

We use periodic boundary conditions in the azimuthal di- 
rection. For the velocity field we apply stress-free boundary 
condition at the radial and latitudinal boundaries. The mag- 
netic field follows a perfect conductor condition at the lower 
radial and the two latitudinal boundaries. On the outer radial 
boundary, we force the field to be radial; for discussion on 
the applicability of this boundary condition for the Sun see 
IWarnecke et al.l (|2Q12b|) . We set the gradient of the temper- 
ature at the lower radial boundary, corresponding to a fixed 
radiative flux, and at the radial outer boundary we set the tem- 
perature to a constant value. At the latitudinal boundaries, we 
impose a vanishing first derivative of entropy to suppress a 
heat flux through the boundary. 

Our runs are characterized by the values of the fluid and 
magnetic Reynolds numbers. Re = Urms/^h and Rm = 
^rms/^^f, respectively, where Urms is the volume averaged 
rms velocity in the convection zone, and = 2t:/{R — 
ro) ~ 21/i? is used as a reference wavenumber. To rep- 
resent the turbulent velocities in a proper way, we define 

^rms 

= Y^3/2(ii^ + ul)e(f)r<R where we neglect the differ- 
ential rotation dominated ^-component. We also define the 
fluid and magnetic Prandtl numbers Pr = i^/xsGS = 
u/t] = Rm/Re, the Coriolis number Co = 2fto/urmskf 
and Taylor number Ta = {2QoR'^ /u)'^. Time is given in 
turnover times, r = {urmsh)~^- We measure the magnetic 
field strength as the rms value averaged over the convection 
zone Brms, and we normalize this value with the equipartition 
value of the magnetic field defined by B'^^ = iiQ{pu^)r<R- 
We use the (semi-) turbulent Rayleigh number Rat from the 



thermally relaxed state of the run. 

To monitor the solutions in the convection zone, we use two 
different heights, one near the surface at ri = 0.97 i?, and 
one in the middle of the convection zone at r2 = 0.84 i?. 
We use the Pencil Code^ with sixth-order centered finite 
differenc es in space and a th ird-order Runge-Kutta scheme in 
time; see iMitra etaH (l2009h for the extension of the Pencil 
Code to spherical coordinates. 

3. RESULTS 

We focus in this work on three simulations which are sum- 
marized in Table [T] The main differences between these 
runs is their rotation rate and the magnetic Reynolds number. 
Run Ac is a continuation of Run A after t/r = 1 150, but with 
a smaller diffusivity rj in the convection zone. The Runs A 
and Ac have a higher Coriolis number Co and lower fluid and 
magnetic Reynolds number Re and Rm than Run B. We show 
the time evolution of the total rms velocity and magnetic field, 

averaged over the whole domain, u\^^ = {u^ -\- Uq -\- 

and 5jms = (^r + + ^')r£' '^^ ^11 the three runs in 
Figure [21 Here, the subscripts on angle brackets denote av- 
eraging over r, 6>, (/). Convection is sufficiently supercritical 
to develop during the first few tens of turnover times. Af- 
ter 50-200 turnover times, the dynamo starts to operate and a 
magnetic field grows at a rate that is higher for faster rotation 
(Run A). Due to high rotation rate and the lower density in the 
corona, the velocities there grow to higher values than in the 
convection zone. As we describe in Section^ we use a higher 
viscosity to suppress these velocities and associated numer- 
ical difficulties. After the magnetic field in the convection 
zone has reached sufficient strength and expanded through- 
out the whole domain, it quenches the high velocities in the 
corona significantly, as evident from Figure |2] When the mag- 
netic field reaches ^J^s ^ 0-35eq the rms velocity decreases 
from iij:^g ^ 6 to ^ 1, i.e., the contribution from the corona 
is now sub-dominant. This is caused by the Lorentz force, 
which becomes much stronger and comparable to the Corio- 
lis force in the corona when the magnetic field has expanded 
into the whole domain. In the saturated state u^^^^ ~ ^rms, 
which is reached after around t/r = 1000 turnover times for 
Runs A and Ac, even though the magnetic field is still grow- 
ing slowly. For Run B, at first it seems that the saturated state 
has been reached at t/r = 1000, but it turns out that both 
^Jms magnetic field start growing again. While the dif- 
ferential rotation profile remains roughly unchanged despite 
of the growth of the energies (see Section [iJl . the magnetic 
field seems to undergo a mode change from oscillatory into 
a stationary solution or a n osc illatory solution with a much 
longer period (see Section [T2l) . 

In Figure [3] we show the balance of various energy fluxes, 
contributing to the total luminosity of the simulated star for 
Run A. The radial components of radiative, convective, ki- 
netic, viscous, and Poynting fluxes, as well as the flux due to 
the turbulent heat conductivity, are defined as 

7'rad = -<i^r^), (12) 
F,on.=Cv{{pUryT'), (13) 

^ http : / /pencil-code . google code . com 
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TABLE 1 
Summary of the runs. 



Run 


grid 


Pr 


Pm 


Ta 


Po/Ps 


Ma 


Rat 


Re 


Rm 


Co 


-^rms/ -^eq 


An 


A 


400 X 256 X 192 


5 


1 


1.4- 10^0 


14 


0.08 


2.6 • 10^ 


25 


25 


11 


0.29 


-0.01 


Ac 


400 X 256 X 192 


5 


1.67 


1.4 • 10^0 


14 


0.08 


2.6 • 10^ 


25 


42 


11 


0.27 


-0.03 


B 


400 X 256 X 192 


4 


1 


7.2 • 10^ 


15 


0.09 


2.6 • 10^ 


36 


36 


5.3 


0.37 


-0.05 



Note. — The second to fourth columns show quantities that are input parameters to the models whereas the quantities in the last four columns are 
results of the simulations computed from the saturated state. The Mach number is defined as Ma = nrms/cs |r=o.97 r and latitudinal differential 
rotation is quantified through A = dQ/dcos'^9\r,<R/Qo. 
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Fig. 2. — Time evolution of the total urms velocity and magnetic field. The 
rms velocity of the whole domain w^ms normalized by Wrms (solid line) 
and is plotted together with the rms magnetic field of the whole domain B^^^ 
normalized by the equipartition field in the convection zone, Beq, (dotted 
line) and multiplied by ten for visualization purposes for Run A (black line), 
Run Ac (blue), and Run B (red). 

F^in = -HiP^ryu'^). (14) 



^Foyn = {EoB(f) — E(f)Be) / flQ^ 



F^,,, = -2u{{puiySir), (15) 
^SGS = -(i^.^''^), (16) 

(17) 

where E = r]fioJ — u x S is the electric field, the primes 
denote fluctuations, and angle brackets averaging over and 
(f) and a time interval over which the turbulence is statisti- 
cally stationary. The convective flux dominates inside the 
convection zone and reaches here much hig her values than 
in our earlier model (* Wamecke et al.ll2Q12bl) . (In our case, 
XsGS ^ 0-02xto, where xto = Urms/^h is an estimate for 
the macrophysical turbulent diffusivity.) In the corona the 
cooling flux keeps the total flux constant. Note the convec- 
tive overshoot into the exterior and the negative radiative flux 
just above the surface (r = R), caused by the higher tempera- 
ture in the corona. The kinetic energy flux has small negatives 
values in the convection zone. The luminosity due to viscosity 
as well as the Poynting flux is too small to be visible. 

3.1. Differential rotation and meridional circulation 

In Figure m we plot the mean rotation profiles (^(r, 0) = 
^0 + U(f)/r sin 6 together with the meridional circulation 
{ur^ue^O) of Runs A, Ac and B. The contours of constant 
rotation are clearly not cylindrical. They show a 'spoke-like' 
structure, i.e. the contours are more radial than cylindrical, 
which is similar to the solar rotation profile obtained by he- 
lioseismology. The equator is rotating faster than the poles, 
which has been seen in previous simulations of .Brun et al.. 



Fig. 3. — The different contributions to the total luminosity (solid pur- 
ple line) are due to radiative diffusion (dashed red line), resolved convection 
(blue dotted), unresolved turbulent convection (black dotted), viscosity (yel- 
low dashed), flux correspond to the cooling function Fcooi (dash-dotted), and 
the Poynting flux (orange dash-dotted) for Run A. The thin solid black lines 
denote the zero level and the total luminosity through the lower boundary, 
respectively. 



(l2QQl . iMieTch et al.1 (l2006h . and lKapvla et"an ([20TQl.[20TTbh 

and resembles the observed rotation of the Sun for our slower 
rotating case (Run B). 

To produce spoke-like rotation contours, the Taylor- 
Proudman theorem has to be overcome by an important con- 
tribution in the evolution equation for the mean azimuthal vor- 
ticity uJ^, which is given by: 



dt 



dz 



VT X Vs] 



(18) 



where d/dz = cosOd/dr — r~^sinO d/dO is the derivative 
along the rotation axis. The first term in Equation ([TSl l is re- 
lated to the curl of the Coriolis force and vanishes for cylin- 
drical ft contours. The second term is the mean baroclinic 
term, which is caused mainly by latitudinal entropy variations. 
We ignore here additional contributions such as meridional 
Reynolds and Maxwell stresses which turn out to be small. In 
Figure in we plot the first and second terms of Equation (fTSl l 
for each of the three runs. These two contributions balance 
each other nearly perfectly. This lets us conclude that these 
two terms give the dominant contribution to the production of 
mean azimuthal vorticity, and that the Taylor-Proudman the- 
orem is broken by the baroclinic term. It is remarkable that 
there is such a large and spatially coherent latitudinal entropy 
gradient, which is crucial to having a significant azimuthal 
baroclinic term, is produced self-consistently in the simula- 
tions. 

Similar results have so far only been obtained in mean-field 



5 






0.0 F , 

0.6 0.7 0.8 0.9 1.0 1.1 
{r/R)sind 



-0.03 




■?':' 


rv ' ' ' 1 


'Vr""!"! ^ 






+ 15° 1 1 






+ 30° 


:':/ ' 




-7° 






-15° 1 1 






-30° 














;,; : 


■ \ \ / ■ " 






^ ^2 





0.9 



Fig. 4. — Differential rotation and meridional circulation zoomed in the northern hemisphere in the convection zone for the three different Runs A, Ac, and B, 
from the top row to the bottom one. Left column: mean rotation profiles Q{r, 0) /flo color coded and as rotation contours. Middle column: mean rotation profiles 
(color coded) overplotted with vectors (white) showing the meridional circulation in terms of the mass flux 'p(ur,ue, 0)- Note the different plotting region in 
the first and middle columns. The white dotted lines indicate the radii r = ri = 0.97 R and r = r2 = 0.84 R. The three black dash-dotted lines represent 
the latitudes 90° — Oi = 30°, 90° — 62 = 15°, and 90° — 0^ = 7° , which are used in the right column. The white-black dashed line indicates the surface 
(r = R). Right column: latitudinal mass flux 'pug / pQUrms plotted over radius r /R for three different latitudes Oi (blue dashed line), O2 (red solid line) and ^3 
(black solid line), po represents the value at the bottom of the convection zone. The dotted and dash-triple-dotted lines represent the corresponding values in the 
southern hemisphere; see the legend for details. The black dashed lines indicate the zero line, the surface (r = R) and the depths, ri and r2- 



simulations ( Brandenb urg et al.ll 19921: iKitchatinov & Riidigei^ 
[1995) by including an anisotropic convective heat conductiv- 
ity, and in convection simulations by prescribing a latitudinal 
entropy gr adient at the lower radial boundary of the convec- 
tion zone (iMiesch et al.ll2QQ6l) . The last column of Figure [5] 



shows the mean latitudinal entropy gradient Vqs for the three 
runs. The spatial distribution of the gradient agrees with the 
baroclinic term as well, because VrT ^ const, so we can 
conclude that the dominant contribution in the baroclinic term 
is due to the product of the latitudinal entropy gradient and the 
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Fig. 5. — Representation of the two dominant terms in the balance of mean 
azimuthal vorticity together with the latitudinal entropy gradient. In the first 
column the mean azimuthal component of the baroclinic term of the total 
entropy gradient and total temperature gradient (VT x Vs)^ is plotted to- 

2 

gether with advection term r sin ^ dQ / dz in the middle column. The values 
are normalized by Qq. The last column shows the mean latitudinal entropy 
gradient RVes/cp. From top to bottom, Run A, Run Ac, and Run B. The 
white/black dashed line indicate the surface (r = R). 



radial temperature gradient, which is more important than the 
radial entropy gradient multiplying the latitudinal temperature 
gradient. 

We also note that, owing to the coronal envelope, differen- 
tial rotation is able to develop a near- surface- shear layer. In 
the three_simulations we can observe a concentration of con- 
tours of Q near the surface, which indicates the presence of 
strong shear; see Figure |4] However, only in Runs A, and Ac, 
the strong shear is below or at the surface, while in Run B, 
it is concentrated at around r = 1.1 R. In our simulation of 
Runs A and Ac the shear layer is more extended than in the 
Sun and penetrates deeper into the convection zone. Further 



Fig. 6. — Off-diagonal component X9r and the radial component Xrr of 
the turbulent heat conductivity tensor normalized by xto = iirms/3fcf and 
calculated from Equations (|22j and f23\ . From left to right. Run A, Run Ac 
and Run B. Note, that the high values at the bottom of the convection zone, 
are due to the vanishing radial entropy gradient. Otherwise the legend is 
similar to Figure |5] 



Studies using higher stratification should prove if this is just 
an artifact of weak stratification, as one might expect. The 
spoke-like rotation profile with the strong shear near the sur- 
face occurs mostly at lower latitudes, i.e. close to the equator. 
At latitudes above ±30° the contours of constant rotation are 
more complex, but show some indication of strong shear close 
to the surface and only in Run B the contours become cylin- 
drical beyond ±60° latitude. The location of the spoke-like 
differential rotation profile coincides with a similarly shaped 
mean latitudinal entropy gradient. The entropy gradient in the 
northern (southern) hemisphere is negative (positive) below 
±30° latitude. In Run B, this region reaches to higher lati- 
tudes than in Runs A and Ac, which leads to radial contours 
of angular velocity at higher latitudes. 

As mentioned above, we expect that latitudinal entropy gra- 
dient to be a consequence of an anisotropic convective diffu- 
sivity tensor. Such anisotropics are caused by the rotational 
influence on the turbulence. In particular, there is a term pro- 
portional to rtoi^oj (Qoi is the ith component of fto, and not 
to be confused with Q{r^ 0)), which gives a symmetric contri- 
bution XrO = X9r proportional to cos sin 0, so it vanishes at 
the poles and at the equator. In the presence of a latitudinal 
entropy gradient, it would lead to an additional contribution 
to the radial convective flux, 

Fr = -Xrr'pTVrS - XrO'pTV eS. (19) 

Since XrO = XOr, and since there is a radial entropy gradient. 
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it would also lead to a contribution in the latitudinal flux, 

Fe = -xerpTVrS - xeepTVes. (20) 

If we ignore the second term proportional to V6>s, we could 
estimate xor by measuring 

Fe = cpp^, (21) 

so 

Xer^ -cp'^/rWrS. (22) 

The result is shown in Figure |6l where we also plot a similar 
estimate of the radial component, 

Xrr ^ -cp<T^/TV^s. (23) 

We normalize both by Xto = ^rms/3/cf and find that the latter 
is about 3 in those units, and the former one is about 1. 

In reality, and as a consequence, we cannot neglect the sec- 
ond term proportion al to V^g. In fact, as show n in the mean 
field calculations of iBrandenburg et aH (|1992|) . it will try to 
balance the first term so as to reduce the latitudinal heat flux 
and thus to produce a latitudinal entropy gradient and thus a 
baroclinic term as we see it. 

Another important result is the solar-like meridional cir- 
culation in the convection zone. In the middle column and 
right column of Figure |4] we plot the meridional circulation 
in terms of the mass flux as vectors of p{ur, uq^ 0) and as ra- 
dial cuts of 'puo through three values of colatitude 9i = 60°, 
O2 = 75°, and O3 = 83°, corresponding to 30°, 15°, 7° 
latitudes. Note here that in the northern (southern) hemi- 
sphere a negative ~pue means poleward (equatorward) flow, 
and a positive one equatorward (poleward) flow. Runs A and 
Ac shows significant solar-like meridional circulation, while 
Run B more diffuse pattern. Looking at Run A, in the north- 
ern hemisphere at lower latitudes (< 20°) just below the sur- 
face r = ri = 0.97 i?, the meridional circulation is pole- 
ward with pue = — 0.007po^rms- Above the surface there is 
a return flow in the equatorward direction. This return flow 
peaks above the surface with a similar flux. The turning point 
~puo = is just below the surface, at around r = 0.985 R. If 
we were to redefine the surface of the simulated star to this 
radius, we would obtain a solar-like meridional circulation. 
In the southern hemisphere the pattern is the same, however 
the meridional circulation is a bit weaker. This behavior can 
be found also in Run Ac, where the flows are weaker and the 
turning point lies slightly deeper. Note that the larger return 
flow is a cause of our particular setup, which has a much lower 
stratification than the Sun, see Figure [T] Higher stratification 
should lead to a much weaker return flow above the surface. In 
Run B, a poleward flow develops in the northern hemisphere 
only close to equator (0 = 6^3) or deeper in the convection 
zone at high latitudes (6 = 61). However, in the southern 
hemisphere, the pattern is similar to the slower rotating runs, 
but with weaker circulation. The meridional circulation has 
a latitudinal dependence. In Runs A and Ac the return flow 
reaches higher velocities at lower latitudes (< ±20°). The 
same is valid for the actual poleward circulation below the 
surface. In Run B, we find the opposite and the return flow in- 
creases with latitude, and so does the meridional circulation. 

In the right panel of Figure lU we can also investigate the 
amount of meridional circulation cells for the lower latitudes. 
We find that there are at least two cells in the convection zone. 
In Run A, there is one cell with poleward flow maxima around 
r = ri and minima around r = 0.91 i?, where the mass flux 



closer to the surface is as large as in the return flow in the 
convection zone. A second one is deeper down in the con- 
vection zone and has a similar extent and flux. In Run Ac, 
the pattern is different and the flux is weaker, but there seems 
to be an indication of a third cell. Note that the stratification 
leads to stronger mass fluxes over a smaller cross-sectional 
area deeper in the convection zone than near the surface, while 
the velocity is similar. In Run B, there are two strong cells 
of meridional circulation. The cells seem to be more cylin- 
drical than latitudinal, which can also be seen in the phase 
shift of the pattern for different latitudes. Even deeper down 
in the convection zone, the meridional flow is much stronger 
than in Runs A and Ac. This is not surprising: rapidly rotat- 
ing stars are expecte d to have a we aker meridional circulation 
than slower rotators (lKohlerill97Qh . 

3.2. Mean magnetic field evolution 

The turbulent helical motions generated by convective heat 
transport together with differential rotation produce a large- 
scale magnetic field inside the convection zone. It grows 
exponentially and shows an initial saturation after around 
t/r = 100 for Run A, where it is still growing but slowly 
and linearly. Run B shows a more peculiar behavior; the field 
seems to have saturated at around t/r = 300, but it starts 
growing again at around t/r = 700 and has not yet satu- 
rated. The second growth is possibly related to the change 
of the oscillatory mode into stable one; see Figure |7] and be- 
low. The magnetic and fluid Reynolds numbers of Run B 
are higher, which should lead to a higher growth rate. How- 
ever, the rotation rate measured by Co is around half that of 
Run A, which leads to a slower amplification of the field. 
At a later time, around t/r = 1000, the field of Run B be- 
comes comparable to or even stronger than that of Run A. The 
Brms reaches around 0.5 of the equipartition field strength in 
Runs A, and Ac, and 0.6 in Run B; see Tableffl The is 
around 20% lower, because the field is mainly concentrated 
in t he convection zone . Comparing these results with those 
of Kapyla et al.l (1201 2l) without a corona but otherwise com- 
parable setup, the magnetic field in the current simulations is 
significantly weaker. Additionally the growth rate of the dy- 
namo is greater than in the models without the corona, where 
it takes up to five time longer to reach dynamically important 
field strengths. This is not surprising because the dynamo in 
the two-layer model is less restricted and has more freedom 
for different dynamo modes to be excited. There is no re- 
striction due to the magnetic boundary at the surface, which 
is open in our simulations and restricted to vertical field in 
the convection zone simulations of Kapyla et al. ( 2012). This 
could explain the fast growth in the beginning, but not the 
decreased saturation level. At the moment the field is weaker, 
but in particular Run B shows a strong growth starting again at 
a later time, which is comparable with the time of the runs of 
Kapyl a et~aD ([2012'). This might lead to a new saturation level 
that is as strong as in the runs of Kapyla et al. (2012). On the 
other hand, the runs in this work and the runs of Kap yla et al.l 
(1201 2l 1201 3i) have also differences in other parameters, such 
as stratification, rotation rate Co or Reynolds numbers Re, 
Rm, so direct comparison might not be always possible. 

Looking just at the ^rms in the convection zone or the B^^^ 
in the whole domain in Figure O we find no clear evidence of 
cyclic behavior of the field. However, investigating the dif- 
ferent components of the mean magnetic field we find signs 
of clear oscillatory behavior. In Figure |7] we plot the az- 
imuthal mean magnetic field 5^ and the radial mean field 
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Fig. 7. — Time evolution of the mean magnetic field in the convection zone. From left to right we show the Runs A, Ac, and B and from top to bottom B^j^ and 
Br 2Lir = ri = 0.97R and at r = r2 = O.MR. Dark blue shades represent negative and light yellow positive values. The dashed horizontal lines show the 
location of the equator at ^ = 7r/2. The magnetic field is normalized by its equipartition value, Beq. 








Fig. 8. — Time series of eight snapshots of the mean azimuthal magnetic field B^j) separated by 22 turnover times and covering one full magnetic cycle. Dark 
blue shades represent negative and light yellow positive values. The dashed line indicates the surface (r = R). The field is normalized by the equipartition field 
strength Beq. 



Br over time and latitude at the different heights ri and r2 
for Runs A, Ac, and B. The structure of the magnetic field 
changes as the dynamo evolves from the kinematic regime, 
where the magnetic field is weak and does not significantly 
influence the flow. In Run A, the azimuthal and radial mean 
fields migrate poleward in the kinematic regime close to the 
equator. The cycle frequency is small, just around t/r = 20. 
In Run B, we found a similar behavior. The fast poleward mi- 
gration happens at lower latitudes (±40°) for both runs. We 
recall that in Run A after a short time (t/r ^ 100) the field is 
strong enough to backreact on the flow. Now two things hap- 
pens simultaneously: an oscillating mean magnetic field is 
starting migrating equatorward at higher latitudes and the fast 
poleward migration becomes slower. The period of the equa- 



torward oscillation is longer and is between t/r = 100 and 
150. The poleward migration near the equator slows down un- 
til it finally turns into an equatorward migration aligned with 
the migration at the higher latitudes (t/r = 500). Thus we 
have equatorward migration of the mean radial and azimuthal 
fields at all latitudes until around t/r = 1500. At this time 
the dynamo mode changes, and consequently its latitudinal 
migration pattern changes. The equatorward migrating and 
oscillating field near the poles stays stable during whole sim- 
ulation, but near the equator the field changes. In the northern 
hemisphere there is a transient poleward migration, which is 
in phase with the equatorward migration near the poles. In the 
southern hemisphere the equatorward migration is still dom- 
inant, but a stationary mode is superposed near the equator. 
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The migration patterns are not just features appearing close 
to the surface, but they penetrate the whole convection zone 
until the bottom, as seen in the two bottom rows of Figure [7] 
Therefore, we can rule out that the meridional circulati on is 
the main driver of this migration. As discussed in Section ITTl 
the meridional circulation shows strong variability in radius 
and has at least two cells. 

Run Ac has been started from a snapshot of Run A after 
around t/r = 1150 by increasing the magnetic Prandtl num- 
ber; see Table [T] The pattern of the mean magnetic field is 
not strongly effected by this change of parameter and is simi- 
lar to Run A. There is a clear equatorward migration near the 
poles and some indication of poleward migration and equator- 
ward migration near the equator in the northern hemisphere. 
The stationary fields superimposed to this cycle are similar to 
those in Run A. 

In Run B, where the fluid and magnetic Reynolds numbers 
are higher and the rotation slower than in Run A, the struc- 
ture of the mean field evolution shows some differences. In 
the kinematic regime, the field is similar to that of Run A in 
which it migrates poleward at lower latitudes. Also, as the 
field gets stronger it begins to migrate from higher latitudes 
toward the equator and the fast poleward branch becomes 
slower. The main difference to Run A is that the poleward 
migration does not turn into equatorward migration near the 
equator. In Runs A and Run Ac, the field strength has no clear 
latitudinal dependence. But in Run B the field strength near 
the poles is around the half the strength than near the equator. 
Only during late times the polar branch increases in strength. 
In Runs A and Ac, the radial and azimuthal components have 
approximately the same strength, whereas in Run B, the ra- 
dial mean magnetic field seems to be weaker by a factor of 
two. Also, Run B shows no clear radial dependence in the 
structure of the mean field. The period of the equatorward 
migrating field is around t/r = 200, which is a bit longer 
than for Run A. The poleward migration near the equator has 
an irregular oscillation and is usually not in phase with the 
equatorward migration near the poles. At around t/r = 1000 
the dynamo mode changes significantly. Not only the mag- 
netic field starts to grow (see Figured, but also the magnetic 
field changes from an oscillatory pattern to a stationary one, 
or at least an oscillatory one with a much longer period; see 
Figure [71 In particular in the norther hemisphere the mean az- 
imuthal field shows are strong increase in strength. The field 
pattern now consists of a strong time-independent field with 
a latitudinal dependence. Near the surface (at r = ri; see 
Section |2]) the field seems to migrate slowly toward the equa- 
tor, but it is not possible to identify a migrating pattern in the 
present run. 

If one translates the cycle period of the equatorward migra- 
tion to solar values using of turnover time r of 1 month, we 
obtain a cycle period of 12 years and 16 years for Runs A and 
B, respectively. This would be a typical value in the middle 
of the convection zone. 

To investigate the equatorward migration in the first half of 
the simulations, we plot the mean azimuthal magnetic field 
B(f) for eight different times for Run A, resolving one cy- 
cle; see Figure [H The field penetrates the whole convection 
zone and has up to four regions with different polarities in one 
hemisphere. These polarities are migrating toward the equa- 
tor. In the northern hemisphere at around 45° latitude, there 
is a magnetic field concentration with positive polarity. After 
At/r = 87 (5th panel), at the same location we find a nega- 



tive magnetic field concentration, and again after At/r = 65 
(8th panel) the same polarity appear as in the beginning of the 
cycle. One can see a clear cyclic equatorward migration of the 
field, but it is irregular. The two hemispheres do not show the 
same magnetic field strength and it seems that, from time to 
time, there is only one dominant polarity in one hemisphere, 
while in the other there are three. Note also the strong nega- 
tive magnetic fields near and above the surface, which seem 
to show also an cycle dependency. 

4. CONCLUSIONS 

We have used a model that combines the turbulent con- 
vective dynamo with a coronal layer to reproduce properties 
of the Sun. We found a solar-like differential rotation with 
roughly radial contours of angular velocity at low latitudes. 
This is accompanied by poleward meridional circulation just 
below the surface and a return flow above the surface. Addi- 
tionally the differential rotation profiles show a near- surface- 
shear layer in two of the three simulations we perform. In the 
third simulation, this layer is still exists but it is now located 
above the surface. We identify as the main driver of spoke-like 
differential rotation the non-zero latitudinal entropy gradient, 
which is self-consistently generated. 

The mean magnet ic field show s a pa t tern o f equatorward 
migration similar to Ka pyla et"an (|2012L |2013|) . This pattern 
is mostly visible at higher latitudes and only in two of the 
simulation to lower latitudes. However, in later stages of the 
simulation the equatorward migration is only visible at high 
latitudes and at lower latitudes poleward migration or station- 
ary modes occur. In one of the simulation the dynamo mode 
changed completely to stationary mode on all latitudes, at 
later stages. The dynam o has a shorter excitation time than 
in the earlier work of Ka pyla et"an (l2012b . 

This work leads to the conclusion that the inclusion of a 
coronal layer to a convective dynamo simulation has an influ- 
ence on the fluid and magnetic properties of the interior. In 
recent simulations, we were able to produce recurrent coro- 
nal ejections from the solar surface (IWamecke et al.ll2Q12b|) 
using a two layer approach. In earlier models of forced turbu- 
lence with a co ronal layer ( Warnecke & Brandenburg 20Tol 
IWamecke et al.ll2Ql IL l2Q12a|) . we also found the ejection of 
magnetic helicity out of the dynamo region. These ejections 
can support and amplify the magnetic field due to enhanced 
helicity fluxes. 

Here, we present evidence that even the fluid properties in 
the bulk of the convection zone might be influenced by the 
coronal layer. Solar-like rotation profiles could not be re- 
produced by earlier direct nume rical simulation of convec- 
tive dynamo (iKap yla et al.l l2012h without prescribing a lati- 
tudinal entrop y gradient at the bottom of the convection zone 
(iMiesch et al.l 12006). However, to have more convincing ev- 
idence in support of this, we need to perform a detailed pa- 
rameter study using different coronal sizes and compare them 
with simulations without a corona. 

Another extension of our work would be the measurement 
of magnetic helicity fluxes through the surface and how their 
dependence on the coronal size. To investigate the mechanism 
of the equatorward migration, which is crucial for understand- 
ing the solar dynamo, one should measure the turbulent trans - 
port coefficients throu gh approaches like the test-field method 
(ISchrinner etani2QQ7h . 

In further work we will investigate the possibility of pro- 
ducing coronal ejections using the setup of these runs. In com- 
parison with Warnecke et al.. (2012h) . we use here a corona 
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with a much higher temperature, and a lower plasma beta. It 
will be interesting to see how the coronal ejections are influ- 
enced by this. 
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